Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:

# Manipulación de estructuras
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# Para realizar plots
library(ggplot2)
library(plotly)

## 
## Adjuntando el paquete: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

# Manipulación de varios plots en una imagen.
library(gridExtra)

## 
## Adjuntando el paquete: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

Elimine eval=FALSE para ejecutar las celdas.

Pregunta 1: Estimadores.

Considere una serie de variables aleatorias \(X_i\) que siguen una distribución de Bernoulli de parámetro \(p=0.5\) y un estimador de \(p\) dado por \(\hat{p}_{n, \sigma} = \epsilon_{\sigma} + \frac{1}{n} \displaystyle{\sum_{i=1}^{n}}X_{i}\) donde \(\epsilon_{\sigma} \sim \mathcal{N}(0,\sigma)\). Para \(\sigma = 0, 1, 2, 4\) grafique el valor de \(\hat{p}_{n, \sigma}\) y comparelo con el valor verdadero. Calcule el sesgo del estimador para cada valor de \(\sigma\).

# Función para generar los estimadores
calcular_estimador <- function(n, sigma) {
  # Simulación de las X_i ~ Bernoulli(0.5)
  X <- rbinom(n, 1, 0.5)
  # Ruido epsilon_sigma ~ N(0, sigma)
  epsilon_sigma <- rnorm(1, mean = 0, sd = sigma)
  # Estimador
  p_gorro <- epsilon_sigma + (1/n) * sum(X)
  return(p_gorro)
}

# sigma = 0
cte <- rep(0.5,1000)
x <- 1:1000

# Dataframe para almacenar los resultados
resultados_0 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())
for (n in x) {
  # Calcular el estimador usando la función
  p_gorro <- calcular_estimador(n, 0)
  
  # Almacenar los resultados
  resultados_0 <- rbind(resultados_0, data.frame(n = n, sigma = 0, p_gorro = p_gorro))
}

# Graficar
ggplot(resultados_0, aes(x = n, y = p_gorro)) +
  geom_line(color = "blue") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimador de p (sigma = 0)",
       x = "n",
       y = expression(hat(p)[n, 0])) +
  theme_minimal()

# Sesgo para sigma = 0
sesgo_estimador_sigma_0 <- mean(resultados_0$p_gorro) - 0.5
print(sesgo_estimador_sigma_0)

## [1] -0.0006492864

# sigma = 1
cte <- rep(0.5,1000)
x <- 1:1000

# Dataframe para almacenar los resultados
resultados_1 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())

for (n in x) {
  # Calcular el estimador usando la función
  p_gorro <- calcular_estimador(n, 1)
  
  # Almacenar los resultados
  resultados_1 <- rbind(resultados_1, data.frame(n = n, sigma = 1, p_gorro = p_gorro))
}

# Graficar
ggplot(resultados_1, aes(x = n, y = p_gorro)) +
  geom_line(color = "green") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimador de p (sigma = 1)",
       x = "n",
       y = expression(hat(p)[n, 1])) +
  theme_minimal()

# Sesgo para sigma = 1
sesgo_estimador_sigma_1 <- mean(resultados_1$p_gorro) - 0.5
print(sesgo_estimador_sigma_1)

## [1] -0.03590337

# sigma = 2
cte <- rep(0.5,1000)
x <- 1:1000

# Dataframe para almacenar los resultados
resultados_2 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())

for (n in x) {
  # Calcular el estimador usando la función
  p_gorro <- calcular_estimador(n, 2)
  
  # Almacenar los resultados
  resultados_2 <- rbind(resultados_2, data.frame(n = n, sigma = 2, p_gorro = p_gorro))
}

# Graficar
ggplot(resultados_2, aes(x = n, y = p_gorro)) +
  geom_line(color = "red") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimador de p (sigma = 2)",
       x = "n",
       y = expression(hat(p)[n, 1])) +
  theme_minimal()

# Sesgo para sigma = 2
sesgo_estimador_sigma_2 <- mean(resultados_2$p_gorro) - 0.5
print(sesgo_estimador_sigma_2)

## [1] -0.02068251

# sigma = 4
cte <- rep(0.5,1000)
x <- 1:1000

# Dataframe para almacenar los resultados
resultados_4 <- data.frame(n = numeric(), sigma = numeric(), p_gorro = numeric())

for (n in x) {
  # Calcular el estimador usando la función
  p_gorro <- calcular_estimador(n, 4)
  
  # Almacenar los resultados
  resultados_4 <- rbind(resultados_4, data.frame(n = n, sigma = 4, p_gorro = p_gorro))
}

# Graficar
ggplot(resultados_4, aes(x = n, y = p_gorro)) +
  geom_line(color = "purple") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimador de p (sigma = 4)",
       x = "n",
       y = expression(hat(p)[n, 1])) +
  theme_minimal()

# Sesgo para sigma = 4
sesgo_estimador_sigma_4 <- mean(resultados_4$p_gorro) - 0.5
print(sesgo_estimador_sigma_4)

## [1] -0.03427481

Respuesta

¿Qué puede decir del sesgo de los estimadores \(\hat{p}_{n,\sigma}\)?¿Son los estimadores con menos sesgo más precisos? Justifique.
- Como el sesgo es la diferencia entre el valor esperado del estimador y el verdadero valor del parámetro, al calcular el sesgo para los diversos valores de \(\sigma\), obtenemos que para \(\sigma=0\) y \(\sigma=1\) el sesgo es cercano a 0. El ruido añadido es muy pequeño o inexistente, por lo que el estimador está estimando correctamente el parámetro \(\hat{p}_{n,\sigma}\).
- Para \(\sigma=2\), el sesgo es negativo, lo que sugiere que el estimador está subestimando el valor de \(\hat{p}_{n,\sigma}\).
- Para \(\sigma=4\), el sesgo es positivo, lo que significa que el estimador está sobreestimando el valor de \(\hat{p}_{n,\sigma}\).
- Por lo tanto los estimadores con menor sesgo tienden a ser más precisos, ya que en promedio sus valores están más cerca del verdadero valor de \(\hat{p}_{n,\sigma}\).

Pregunta 2: Intervalo de Confidencia

En las preguntas 2, 3 y 4 deberá trabajar con el dataset diabetes_prediction_dataset.

datos <- read.csv("diabetes_prediction_dataset.csv", header = TRUE)
head(datos)

Una forma de estimar la distribución de la media de una población es a través de la realización de la sampling distribution de una distribución cualquiera. El valor obtenido en cada sampling distribution nos entrega un estimador de la media que posee una determinada distribución de la población. Sabiendo esto, calcule la sampling distribution de las columnas bmi y blood_glucose_level, obteniendo el intervalo de confianza de \(95\%\) para cada una de las medias obtenidas desde la distribución. Para realizar este experimento considere los siguientes puntos:

Obtener las medias y desviaciones estándar de la población real.
Realizar una sampling distribution con un tamaño de muestra igual a \(100\). Repita la obtención de la media un número elevado de veces (recomendación \(5000\) veces).
Calcular el intervalo de confianza para cada uno de las medias obtenidas en las iteraciones.
De acuerdo a los valores obtenidos (medias e intervalos de confianza), grafique cada una de las medias obtenidas en conjunto a sus intervalos de confianza. Aquí debe notar que, si el intervalo de confianza contiene a la media de la población, este se considerará como parte del intervalo de confianza del \(95\%\), haga un conteo de cuantos valores están contenidos en este intervalo. Para graficar los intervalos deconfianza puede ser útil utilizar forest plot.
Compare las estimaciones para bmi y blood_glucose_level y señale que diferencias y similitudes encuentra entre estas (en cuanto a la calidad de la estimación, no los valores obtenidos). De una explicación de por qué cree que se dan las similitudes/diferencias.

Hints:

Para comparar comparar las 2 estimaciones puede ser útil normalizar bmi y blood_glucose_level a una misma escala.
Para realizar la sampling distribution podría serle útil el comando sample().
Puede ser útil generar un dataframe para graficar con ggplot2.

Respuesta:

# Definimos tamaños de muestreo
tamano_muestra = 100
n_muestras = 5000

# Inicializamos estructuras necesarias
list_mean_bmi = vector()
list_interval_bmi = vector()
list_typeCI_bmi = vector()
list_prop_bmi = vector()
sucess_bmi = 0

list_mean_bgl = vector()
list_interval_bgl = vector()
list_typeCI_bgl = vector()
list_prop_bgl = vector()
sucess_bgl = 0

#función para calcular desviación estándar
sd.p <- function(x){sd(x)*sqrt((length(x)-1)/length(x))}

# Obtenemos la media y desviación estándar de cada columna
mu_bmi <- mean(datos$bmi)
sd_bmi <- sd.p(datos$bmi)
mu_bmi

## [1] 27.32077

sd_bmi

## [1] 6.63675

mu_glucose <- mean(datos$blood_glucose_level)
sd_glucose <- sd.p(datos$blood_glucose_level)
mu_glucose

## [1] 138.0581

sd_glucose

## [1] 40.70793

alpha <- 0.05

# Sampling distribution, calculo del intervalo de confianza y proporción.

for (i in 1:n_muestras) {
  ### Para bmi ###
  samples_bmi <- sample(datos$bmi, tamano_muestra, replace = TRUE)
  sample_mu_bmi <- mean(samples_bmi)
  sample_sd_bmi <- sd(samples_bmi)
  
  # Calculo de intervalo de confianza
  se_bmi <- sample_sd_bmi/sqrt(tamano_muestra)
  error_bmi <- qnorm(1-alpha/2) * se_bmi
  left_bmi <- sample_mu_bmi - error_bmi
  right_bmi <- sample_mu_bmi + error_bmi
  
  # Se guardan los resultados 
  list_mean_bmi[i] <- sample_mu_bmi
  list_interval_bmi <- c(list_interval_bmi, left_bmi, right_bmi)
  list_prop_bmi[i] <- sucess_bmi/i
  
  # Se verifica si la media de la población se encuentra dentro del intérvalo de confianza
  if (left_bmi <= mu_bmi && right_bmi >= mu_bmi) {
    sucess_bmi <- sucess_bmi + 1
  }
  
  
  ### Para blood_glucose_level ###
  samples_glucose <- sample(datos$blood_glucose_level, tamano_muestra, replace = TRUE)
  sample_mu_glucose <- mean(samples_glucose)
  sample_sd_glucose <- sd(samples_glucose)
  
  # Calculo de intervalo de confianza
  se_glucose <- sample_sd_glucose / sqrt(tamano_muestra)
  error_glucose <- qnorm(1 - alpha / 2) * se_glucose
  left_glucose <- sample_mu_glucose - error_glucose
  right_glucose <- sample_mu_glucose + error_glucose
  
  # Se guardan los resultados
  list_mean_bgl[i] <- sample_mu_glucose
  list_interval_bgl <- c(list_interval_bgl, left_glucose, right_glucose)
  list_prop_bgl[i] <- sucess_bgl / i
  
  # Se verifica si la media de la población se encuentra dentro del intérvalo de confianza
  if (left_glucose <= mu_glucose && right_glucose >= mu_glucose) {
    sucess_bgl <- sucess_bgl + 1
  }
  
}

# Generamos dataframes para plotear mas facilmente los datos.

# Convertimos los vectores de intervalos en matrices de dos columnas
list_interval_bmi_matrix <- matrix(list_interval_bmi, ncol = 2, byrow = TRUE)
list_interval_bgl_matrix <- matrix(list_interval_bgl, ncol = 2, byrow = TRUE)

# Generamos dataframes para graficar más fácilmente los datos
df_bmi <- data.frame(media = list_mean_bmi, ci_low = list_interval_bmi_matrix[, 1], ci_high = list_interval_bmi_matrix[, 2], prop_ci = list_prop_bmi)
df_bgl <- data.frame(media = list_mean_bgl, ci_low = list_interval_bgl_matrix[, 1], ci_high = list_interval_bgl_matrix[, 2], prop_ci = list_prop_bgl)

# Verificar si el intervalo de confianza contiene la media de la población
df_bmi$contiene_media <- (df_bmi$ci_low <= mu_bmi & df_bmi$ci_high >= mu_bmi)
df_bgl$contiene_media <- (df_bgl$ci_low <= mu_glucose & df_bgl$ci_high >= mu_glucose)

# Plot de Intervalos de confianza
# Graficar una muestra aleatoria de 100 medias e intervalos de confianza para BMI
set.seed(123)  # Para reproducibilidad
sample_indices_bmi <- sample(1:n_muestras, 100)  # Seleccionar 100 índices aleatorios
sample_bmi <- df_bmi[sample_indices_bmi, ]

# Crear gráfico para la muestra de BMI
fp_bmi <- ggplot(data = sample_bmi, aes(x = factor(1:100), y = media, ymin = ci_low, ymax = ci_high)) +
  geom_pointrange() + 
  geom_hline(yintercept = mu_bmi, lty = 2) +  # Línea horizontal para la media poblacional
  coord_flip() +  # Girar coordenadas
  xlab("Muestra Aleatoria de Medias BMI") + ylab("Media (95% CI)") +
  theme_bw() + 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())  # Opcional: limpiar etiquetas del eje y
print(fp_bmi)

# Plot de Intervalos de confianza
# Graficar una muestra aleatoria de 100 medias e intervalos de confianza para Glucose
set.seed(123)  # Para reproducibilidad
sample_indices_bgl <- sample(1:n_muestras, 100)  # Seleccionar 100 índices aleatorios
sample_bgl <- df_bgl[sample_indices_bgl, ]

# Crear gráfico para la muestra de Glucose
fp_bgl <- ggplot(data = sample_bgl, aes(x = factor(1:100), y = media, ymin = ci_low, ymax = ci_high)) +
  geom_pointrange() + 
  geom_hline(yintercept = mu_glucose, lty = 2) +  # Línea horizontal para la media poblacional
  coord_flip() +  # Girar coordenadas
  xlab("Muestra Aleatoria de Medias Glucose") + ylab("Media (95% CI)") +
  theme_bw() + 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())  # Opcional: limpiar etiquetas del eje y
print(fp_bgl)

# Suponiendo que ya tenemos los datos de las muestras y las estadísticas
# Normalizamos las variables
normalizar <- function(x) {
  return((x - mean(x)) / sd(x))
}

# Normalizamos las medias obtenidas de las muestras
list_mean_bmi_normalized <- normalizar(list_mean_bmi)
list_mean_bgl_normalized <- normalizar(list_mean_bgl)

# Creamos un data frame para facilitar el plotting
df_comparacion <- data.frame(
  index = 1:n_muestras,
  bmi = list_mean_bmi_normalized,
  blood_glucose_level = list_mean_bgl_normalized
)

# Convertimos el data frame a formato largo sin reshape2
df_long <- data.frame(
  index = rep(df_comparacion$index, 2),
  Variable = rep(c("BMI", "Blood Glucose Level"), each = n_muestras),
  Valor = c(df_comparacion$bmi, df_comparacion$blood_glucose_level)
)

# Graficamos usando ggplot2
library(ggplot2)

ggplot(df_long, aes(x = index, y = Valor, color = Variable)) +
  geom_line() +
  labs(title = "Comparación de Estimaciones Normalizadas de BMI y Glucosa en Sangre",
       x = "Índice de Muestra",
       y = "Estimación Normalizada") +
  theme_minimal() +
  theme(legend.title = element_blank())

# Plot de proporción de CI
p_prop_bmi <- ggplot(df_bmi, aes(x = 1:n_muestras, y = prop_ci)) +
  geom_line(color = "blue") +
  ggtitle("Proporción de Intervalos de Confianza que Contienen la Media - BMI") +
  xlab("Iteraciones") + ylab("Proporción Acumulada") +
  ylim(0, 1)

p_prop_bgl <- ggplot(df_bgl, aes(x = 1:n_muestras, y = prop_ci)) +
  geom_line(color = "green") +
  ggtitle("Proporción de Intervalos de Confianza que Contienen la Media - Glucose") +
  xlab("Iteraciones") + ylab("Proporción Acumulada") +
  ylim(0, 1)

# Mostrar las gráficas
print(p_prop_bmi)

print(p_prop_bgl)

Respuesta Es posible observar que la estimación para blood_glucose_level tiene mejores resultados ya que a simple vista se observan más valores cercanos a la media real, en cambio para bmi se observan valores mas dispersos y no tan cercanos a la media real por lo que se espera eque este estimador tenga un mayor error.

Pregunta 3: Estimación de Máxima Verosimilitud

El objetivo de esta parte será obtener los parámetros de las distribuciones de la media y la mediana de bmi. Para responder la pregunta realice los siguientes puntos:

Simule la distribución empírica a través de mustras con repetición, como las de la pregunta anterior, y grafique los histogramas de la media y mediana.
Utilice la estimación logaritmica de máxima verosimilitud negativa para obtener los parámetros de estas distribuciones.
Gráfique a traves de un gráfico de calor el rango de valores en que se mueve la solución del problema de likelihood (puede ser útil persp() o filled.contour()).
Aplique el comando nlminb() sobre la likelihood y encuentre el máximo o mínimo del problema a optimizar.
Genere muestras lo suficientemente grandes de las distribuciones normales con los parámetros obtenidos por máxima verosimilitud, y comparelas con la distribución empírica mediante histogramas.

Cabe señalar que el método de máxima verosimilitud deberá ser programado por usted y no podrá utilizar librerías que entreguen el valor directo (por ejemplo la librería MASS).

Respuesta

# Obtenemos métricas de muestreo con repetición
set.seed(99)
means <- c()
medis <- c()
n_muestras <- 10000 
tamano_muestra <- 100 

for (i in 1:n_muestras) {
  muestra <- sample(datos$bmi, tamano_muestra, replace = TRUE)
  means <- c(means, mean(muestra))
  medis <- c(medis, median(muestra))
}

# dataframe con las medias y medianas de las muestras
resultados <- data.frame(media = means, mediana = medis)

# Histograma media
ggplot(resultados, aes(x = media)) + 
  geom_histogram(binwidth = 0.1, fill = "blue", color = "white") + 
  ggtitle("Histograma de Medias Simuladas") + 
  theme_minimal()

# Histograma mediana
ggplot(resultados, aes(x = mediana)) + 
  geom_histogram(binwidth = 0.1, fill = "green", color = "white") + 
  ggtitle("Histograma de Medianas Simuladas") + 
  theme_minimal()

# Media
# función log likelihood
ll_plot <- function(mu, sigma) {
  likelihood <- -sum(dnorm(means, mean = mu, sd = sigma, log = TRUE))
  return(likelihood)
}

ll_plot <- Vectorize(ll_plot)

# Definir el espacio donde se va a evaluar ll_plot
mu_m <- seq(26, 29, length.out = 10)    # Rango de valores para la media
sigma_m <- seq(0.01, 1, length.out = 10)   # Rango de valores para la desviación estándar

likelihood_m <- outer(X = mu_m, Y = sigma_m, ll_plot)

filled.contour(x = mu_m, y = sigma_m, z = likelihood_m, xlab = expression(mu), 
               ylab = expression(sigma), color = topo.colors, 
               main = "Log Likelihood para la Media")

# Mediana
ll_plot_mediana <- function(mu, sigma) {
  likelihood <- -sum(dnorm(medis, mean = mu, sd = sigma, log = TRUE))
  return(likelihood)
}
ll_plot_mediana <- Vectorize(ll_plot_mediana)
mu_d <- seq(27, 28, length.out = 10)
sigma_d <- seq(0.01, 1, length.out = 10)

likelihood_d <- outer(X = mu_d, Y = sigma_d, ll_plot_mediana)

filled.contour(x = mu_d, y = sigma_d, z = likelihood_d, xlab = expression(mu), 
               ylab = expression(sigma), color = topo.colors, 
               main = "Log Likelihood para la Mediana")

likelihood_mean <- function(param) {
  # Definimos los parametros de entrada de la funcion
  mu <- param[1]
  sigma <- param[2]
  
  # Definimos la likelihood como la suma logarítmica de la función de densidad
  log_likelihood <- -sum(dnorm(means, mean = mu, sd = sigma, log = TRUE))
  
  return(log_likelihood)
}

# Agrrgue el rango donde operará nlminb
result_mean <- nlminb(objective = likelihood_mean, 
                      start = c(mu = 27, sigma = 1),  # Valores iniciales
                      lower = c(mu = 26, sigma = 0),  # Límite inferior
                      upper = c(mu = 28, sigma = 2))   # Límite superior
# Mostrar los resultados
cat("Estimación por Máxima Verosimilitud para la media:\n")

## Estimación por Máxima Verosimilitud para la media:

cat("Mu óptimo:", result_mean$par[1], "\n")

## Mu óptimo: 27.32434

cat("Sigma óptimo:", result_mean$par[2], "\n")

## Sigma óptimo: 0.6633704

likelihood_med <- function(param) {
  # Definimos los parametros de entrada de la funcion
  mu <- param[1]
  sigma <- param[2]
  
  # Definimos la likelihood como la suma logarítmica de la función de densidad
  log_likelihood <- -sum(dnorm(medis, mean = mu, sd = sigma, log = TRUE))
  
  return(log_likelihood)
}

# Agrrgue el rango donde operará nlminb
result_med <- nlminb(objective = likelihood_med, 
                      start = c(mu = 27, sigma = 1),  # Valores iniciales
                      lower = c(mu = 26, sigma = 0),  # Límite inferior
                      upper = c(mu = 28, sigma = 2))   # Límite superior
cat("Estimación por Máxima Verosimilitud para la mediana:\n")

## Estimación por Máxima Verosimilitud para la mediana:

cat("Mu óptimo:", result_med$par[1], "\n")

## Mu óptimo: 27.30507

cat("Sigma óptimo:", result_med$par[2], "\n")

## Sigma óptimo: 0.1025047

# Muestra de medias
mu_media <- 27.32434
sigma_media <- 0.663372
n_muestra <- 10000
muestras_media <- rnorm(n_muestra, mean = mu_media, sd = sigma_media)
df_muestras_media <- data.frame(muestras_media)
# Histogrmas
hist_media_simulada <- ggplot(df_muestras_media, aes(x = muestras_media)) +
  geom_histogram(binwidth = 0.1, fill = "blue", color = "white") +
  ggtitle("Histograma de Medias Simuladas") +
  theme_minimal()
hist_media_simulada

# Muestra de medianas
mu_mediana <- 27.30507
sigma_mediana <- 0.1025046
n_muestra <- 10000
muestras_mediana <- rnorm(n_muestra, mean = mu_mediana, sd = sigma_mediana)
df_muestras_mediana <- data.frame(muestras_mediana)
# Histogrmas
hist_mediana_simulada <- ggplot(df_muestras_mediana, aes(x = muestras_mediana)) +
  geom_histogram(binwidth = 0.1, fill = "green", color = "white") +
  ggtitle("Histograma de Medianas Simuladas") +
  theme_minimal()
hist_mediana_simulada

Al analizar las distribuciones de medias y medianas obtenidas por máxima verosimilitud, se observa lo siguiente:

Medias: La estimación de máxima verosimilitud captura adecuadamente la variabilidad observada en los datos empíricos. Esto se refleja al comparar los histograma de las medias empíricas y el de las medias simuladas, lo que indica que la distribución normal estimada es un buen ajuste para los datos.
Medianas: En el caso de las medianas, no se observa el mismo comportamiento. El histograma de las medianas empíricas muestra una fuerte concentración en un solo valor, lo que indica una baja variabilidad en los datos originales. Esto difiere con las medianas simuladas, que presentan mayor variabilidad, propia de una distribución más normal.

La diferencia en la variabilidad entre las medianas empíricas y las simuladas sugiere que la distribución de los datos empíricos tiene características que no son bien capturadas por la estimación normal, posiblemente debido a una concentración excesiva de valores en un punto específico.

Pregunta 4: Bootstrap I

En esta pregunta se error e inetrvalo de confianza para la varianza de la columna HbA1c_level.

Las actividades por realizar son las siguientes:

Programar el método Bootstrap para calcular el error estándar de la varianza. No se permite la utilización de librerías de bootstrap para esta parte.
Visualizar a través de un de puntos las distintas varianzas obtenidas al realizar el muestreo con Bootstrap y comparar con el valor real.
Visualizar a través de un gráfico el histograma obtenido al realizar el muestreo con Bootstrap y comparar con el valor real.
Obtener el 95% de intervalo de confianza de la estimación.
¿Qué puede se inferir a pertir de los resutados de Bootstrap?

Respuesta:

# Bootstrap
B <- 5000
largo <- length(datos$HbA1c_level)
output <- c()

for (it in 1:B) {
  sample_data <- sample(datos$HbA1c_level, largo, replace = TRUE)
  output[it] <- var(sample_data)
}

# Valor real de la varianza
real_variance <- var(datos$HbA1c_level)
real_variance

## [1] 1.146339

# Gráfico de puntos
plot(output, main = "Varianzas obtenidas con Bootstrap",
     ylab = "Varianza", xlab = "Muestra de Bootstrap")
abline(h = real_variance, col = "red", lty = 2) # Línea roja para la varianza real
legend("topright", legend = "Varianza real", col = "red", lty = 2)

# Histogrma
hist(output, main = "Histograma de varianzas de Bootstrap",
     xlab = "Varianza", breaks = 30, col = "lightblue")
abline(v = real_variance, col = "red", lty = 2) # Línea roja para la varianza real
legend("topright", legend = "Varianza real", col = "red", lty = 2)

# Obtenemos error e intervalos de confianze
se_vars <- sd(output)
z_a2 <- qnorm(1 - alpha/2)
alpha <- 0.05

# límite inferios
l.CI <- quantile(output, alpha/2)
# límite superior
u.CI <- quantile(output, 1-alpha/2)

sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)

## [1] "El intervalo de confidencia de 95% de las varianzas es: (1.136,1.157)"

sprintf('El SE de la varianzas es: (%.3f)', se_vars)

## [1] "El SE de la varianzas es: (0.005)"

Respuesta:

Para el gráfico de puntos para la varianza obtenida con boostrap se puede ver que los puntos están mayormente concentrados en la línea punteada roja, la cual corresponde a la varianza real. Además se nota que hay puntos outliers presentes pero en menor cantidad.

Para el histrograma de la varianza se cómo claramante las varianzas obtenidas con Bootstrap tienen una distribución cercana a la normal centrada en el valor de la varianza real.

Se puede inferir que al aplicar el método bootstrap para error estándar de la varianza, estos resultados tendrán una distribución cercana a la normal centrada en el valor de la varianza real, más aún en el histograma ya se puede apreciar cuál será un valor aproximado de el 95% de intervalo de confianza de la estimación.

Pregunta 5: Bootstrap II

El siguiente código genera una regresión lineal de bmi con respecto a age, es decir, una función lineal de age, \(y(age) = b + m\cdot age\), que pretende determinar el valor de bmi.

linearMod <- lm(bmi ~ age, data=datos)
print(linearMod)

## 
## Call:
## lm(formula = bmi ~ age, data = datos)
## 
## Coefficients:
## (Intercept)          age  
##    23.15536      0.09945

m <- linearMod$coefficients["age"]
b <- linearMod$coefficients["(Intercept)"]

ggplot() +
  geom_point(aes(x=datos$age, y=datos$bmi)) +
  geom_line(aes(x=datos$age, y=datos$age*m+b), color="red") +
  ggtitle("Regresión lineal") +
  ylab("bmi") + 
  xlab("Age")

Repita el proceso de la pregunta 4 para estimar el error e intervalos de confianza de los parámetros de la regresión (\(m\) y \(b\)). ¿Qué indican los resultados de Bootstrap?¿Un valor bajo en el error significa que la regresión es buena?

# Bootstrap
B <- 500
largo <- length(datos$bmi)
output_m <- numeric(B)
output_b <- numeric(B)

for (it in 1:B) {
  # Tomar una muestra con reemplazo de los datos
  sample_indices <- sample(1:largo, largo, replace = TRUE)
  sample_data <- datos[sample_indices, ]
  
  # Ajustar la regresión lineal con la muestra
  bootstrap_model <- lm(bmi ~ age, data = sample_data)
  
  # Guardar los coeficientes en las salidas correspondientes
  output_m[it] <- bootstrap_model$coefficients["age"]
  output_b[it] <- bootstrap_model$coefficients["(Intercept)"]
}

# Mostrar algunos resultados del bootstrap
head(output_m)

## [1] 0.09902192 0.09975384 0.10031168 0.10003177 0.10001440 0.09812545

head(output_b)

## [1] 23.16575 23.11355 23.10619 23.13401 23.11621 23.19340

# Obtenemos error e intervalos de confianze
se_m <- sd(output_m)
z_a2 <- qnorm(1 - alpha / 2)
alpha <- 0.05

# límite inferios
l.CI <- mean(output_m) - z_a2 * se_m
# límite superior
u.CI <- mean(output_m) + z_a2 * se_m

sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)

## [1] "El intervalo de confidencia de 95% de las varianzas es: (0.098,0.101)"

sprintf('El SE de la varianzas es: (%.3f)', se_m)

## [1] "El SE de la varianzas es: (0.001)"

# Obtenemos error e intervalos de confianze
se_b <- sd(output_m)
z_a2 <- qnorm(1 - alpha / 2)
alpha <- 0.05

# límite inferios
l.CI <- mean(output_b) - z_a2 * se_b
# límite superior
u.CI <- mean(output_b) + z_a2 * se_b

sprintf('El intervalo de confidencia de 95%% de los sesgos es: (%.3f,%.3f)', l.CI , u.CI)

## [1] "El intervalo de confidencia de 95% de los sesgos es: (23.156,23.159)"

sprintf('El SE del sesgo es: (%.3f)', se_b)

## [1] "El SE del sesgo es: (0.001)"

Respuesta

Error estándar bajo: Indica que los coeficientes son consistentes a través de las muestras, lo que sugiere un modelo de regresión estable.

Un error bajo no garantiza que el modelo sea bueno en términos predictivos, pero lo hace fiable de que observada relaciones en los datos actuales.

A work by CC6104

Tarea 2

Tarea 2: Frequentist Inference

Integrantes :

Cuerpo Docente:

Fecha límite de entrega: 08/10/2024

Índice:

Objetivo

Instrucciones:

Referencias: